A miRNAs catalogue from third-stage larvae and extracellular vesicles of Anisakis pegreffii provides new clues for host-parasite interplay

Anisakids are widespread marine parasites of medical, veterinary and economic relevance. They infect marine natural hosts but humans can accidentally acquire the fish-borne zoonosis anisakiasis by ingesting infected raw fishes or mollusks. Among the several species described, Anisakis pegreffii is one of the main etiological agent of the disease, in particular in the Mediterranean area. Despite the growing evidence of miRNAs involvement in host-parasite interplay, and the emerging role of exosomal microvesicles in shuttling them between different cell types (and sometime across species), no information on miRNAs from any Anisakis species is presently available. In this study we isolated extracellular vesicles (EVs) released by Anisakis pegreffii infective third-stage larvae (L3) and analyzed by RNA-seq small RNAs from both L3 and EVs. We showed by nanoparticle tracking analysis that L3 release in culture medium particles of size compatible with the one of extracellular vesicles. A catalogue of 156 miRNAs from A. pegreffii was compiled by sequence comparison to evolutionary close species and miRNA prediction software. Using differential expression analysis, we identified a small number of highly abundant miRNAs in larvae and extracellular vesicles fractions whose potential biological relevance may deserve future investigation. Finally, A. pegreffii miRNAs were compared to those described in other parasitic helminths and predicted targets among human genes were searched, suggesting their potential involvement during infection.

. After quality filtering, adapter trimming and size selection (≥ 14 nt), approximately 140 million reads were retained and mapped to the A. simplex genome (AS14). Reads mapping to AS14 (a total of ~ 55 corresponding to the 38.74% of the reads), excluding those representing ribosomal RNAs, were analysed for their size distribution (Fig. 2). In third stage larvae most reads were in the range of 20-24 nt in length, with two peaks at 21 nt and 23 nt. Although these lengths are fully compatible with the mean size of miRNAs, this bimodal distribution is not typical in small RNA-seq studies. For this reason we analyzed the ten most abundant sequences in the 21 nt and 23 nt peaks and found that they all represented miRNAs and accounted for 84% and 78% of the peak content, respectively. Among highly represented miRNAs were miR-100, miR-71, miR-9, miR-5358a and miR-5360. Extracellular vesicles enriched fractions showed a substantially different size distribution with no evident peaks and a higher frequency of reads 18-23 nt in length. Analysis of the most abundant sequences in this range also showed that Table1. RNA-seq results from the sequencing of the small-RNA fraction isolated from infective third-stage larvae of Anisakis pegreffii and its released extracellular vesicles. Data on raw reads, reads passed cutadapt, reads mapping to the Anisakis simplex genome and to rRNAs are reported. Numbers indicate million reads.  www.nature.com/scientificreports/ miRNAs were highly represented, confirming the overall suitability of our RNA-seq libraries for the planned miRNA analysis.

A. pegreffii miRNAs prediction and counting. Given the absence of information on miRNAs from any
Anisakis species, we used two complementary approaches to obtain a list of putative miRNAs from A. pegreffii. First, we retrieved from miRBase 23 all the available miRNAs from the pig roundworm Ascaris suum, which is the species with experimentally validated miRNAs evolutionary closer to Anisakis species. We then used these miRNAs to search for putative orthologues in the Anisakis simplex genome and, using this in silico approach, we obtained a list of 97 hypothetical A. simplex miRNA precursors and 115 mature miRNAs (dataset 1). In a second approach, as described in detail in the method section, reads from our samples were used to search the A. simplex genome using the miRNA prediction software miRDeep*. This way we obtained a list of 150 putative mature (and corresponding precursors) miRNAs from A. simplex (dataset 2). By merging these two datasets, we obtained a final non-redundant list of 206 putative Anisakis mature miRNAs. Alignments of reads from our L3 and EVs samples to this list, considering only miRNAs with reads in two replicates of at least one sample, provided evidence for the expression of 156 mature miRNAs. Among these, 67 (43%, from 46 miRNA precursors) were putative orthologues of A. suum and were renamed according to the current miRNAs nomenclature, as ape-miR-followed by the identification codes used for A. suum. The remaining putative novel 89 miRNAs (57%) were named as novelMiR followed by the identification number assigned by the miRDeep* software. Most of the miRNAs were observed in the L3 samples (126)  Four miRNAs (ape-miR-100a-5p, ape-miR-1-3p, ape-miR-71-5p and ape-miR-9-5p) were by far the most abundant in L3, with > 100,000 mean CPM. The remaining miRNAs could be classified in four additional categories according to CPM values: (i) > 10,000 n = 8; (ii) > 1,000 n = 15; (iii) > 100 n = 33; (iv) ≥ 10 n = 63). Interestingly, ape-miR-100a-5p was also the most abundant miRNAs in EVs (629,841 mean CPM, 8720 mean counts). The 20 most abundant miRNAs identified in L3 are listed in Table 2. The predicted secondary structures for the first three most abundant miRNAs are shown in Fig. 3.
The mature miRNAs correlation and cluster analyses confirmed the overall good quality of replicates, showing a higher homogeneity in the larval sample than in extracellular vesicles sample, which maintain their clustering trend despite the higher variation, probably due to sample normalization effect on triplicates ( Supplementary  Fig. 2). The expression heatmap from matures miRNAs highlighted groups with specific profile signatures, of which a few corresponded to enriched miRNAs in EVs samples (Fig. 4).
Sample-specific miRNA enrichment was evaluated by pairwise comparisons between the two samples (L3 vs EVs): fold change (FC) and false discovery rates (FDR) were calculated to provide statistical validation (Supplementary Table 2). Using as threshold parameters |log2(FC)|> 1 and FDR < 0.05, we found that 38 miRNAs were www.nature.com/scientificreports/ differentially expressed, with 26 upregulated in L3 and 12 in EVs as shown in the Volcano plot (Fig. 5). Number of differentially expressed miRNAs applying progressively more stringent statistical thresholds are shown in Table 3.

Stem and Loop RT-PCR validation of miRNAs.
With the aim to confirm the presence of miRNAs in the samples studied, a list of ten miRNAs has been selected for their experimental validation in Stem and Loop RT-PCR. The list has been elaborated based on abundance category for L3 (three representative of the first categories of abundance, as > 100,000 mean CPM; > 10,000 mean CPM and > 1,000 mean CPM), on the novelty and on their differential expression (upregulated in EVs). All the ten selected miRNAs were successfully amplified and confirmed using Real Time Stem and Loop PCR (novel-miR-19, miR-7-5p, novel-miR-65, novel-miR-27, miR-72-5p, novel-miR-184, miR-100a-1-5p, lin-4-5p, miR-1-3p, novel-miR-131).

Target analysis and seed conservation in parasitic helminths.
The complementarity between the 3′-UTR of the mRNA target and the miRNA seed region (nucleotides 2-8) is a crucial feature of miRNA-mRNA interaction. Identity of the entire miRNA or of the seed region may be indicative of evolutionary conservation of its function.
Comparison of A. pegreffii miRNAs to those from other helminths as Nematoda (A. suum, Trichuris muris, B. malayi, Nippostrongylus brasiliensis, Heligmosomoides polygyrus, Haemonchus contortus), Trematoda (Fasciola hepatica, Schistosoma mansoni, Schistosoma japonicum) and Cestoda (Dibothriocephalus dendriticus, Mesocestoides corti, Taenia crassiceps, Taenia asiatica, Echinococcus multilocularis) showed different levels of conservation. A complete seed conservation was reported between miRNAs from A. pegreffii with those of parasitic nematodes A. suum as shown in Table 4 (85% of the most abundant in L3 and 46% of EVs enriched) followed by B. malayi (75% of the most abundant in L3 and 31% of EVs enriched). Additionally, complete seed conservation and miR-NAs homology were observed with other helminths as trematodes and cestodes as well as with human miRNAs, suggesting a potential functional role conserved across phylogenetically distant taxa.
The most abundant miRNA in both larvae and extracellular vesicles is ape-miR-100a-5p that shows complete homology to miR-100 from other parasitic helminths and humans. It belongs to the miR-10 family, which besides miR-100 also includes miR-51, miR-57 and miR-99 (Rfam RF00104). The predicted putative gene target www.nature.com/scientificreports/ of ape-miR-100a-5p is TRIB2, a crucial gene for regulation of apoptosis and thymocyte cellular proliferation 24 , and its dysregulation was found associated with tumors, including colorectal cancer 25 . Another abundant miRNA is lin-4-5p: originally identified in C. elegans in relation to developmental timing 26 , it shows complete seed identity with miR-125 of several parasitic helminths and humans. The dendrogram obtained from the comparison of the available orthologues showed similarities between ape-lin-4 in A. suum and B. malayi (Fig. 6). miR-125 seems to be involved in fundamental aspects of parasitic infection: F. hepatica miR-125b (fhe-miR-125b) was the most abundant miRNA trafficked by EVs observed in peritoneal macrophages during the early phase of infection. Given the homology with the mammalian miRNA hsa-miR-125b, it has been suggested the hijacking of the miRNA machinery as a parasitic strategy to control host innate cell function 27 . Ape-lin-4-5p showed several interesting potential gene targets: ARID3B involved in crucial cellular processes as transcriptional regulation 28 , STARD13 involved in cell proliferation and tumor suppression and in particular, miR-125b induces metastasis by targeting STARD13 in in-vitro breast cancer cells 29 . Additional gene targets are BMF, related to apoptosis activation 30 and FREM1, a gene related to IL-1 inflammatory response 31 .
Among miRNAs selectively packaged into extracellular vesicles, novel-miR-19 showed no match with other helminths but a match with human miRNAs was observed (Table 4): its putative gene target is FAM83C, which is involved in regulation of MAPK signaling in cancer cells 32 .

Discussion
Investigations based on high-throughput sequencing techniques as genomic, transcriptomic and proteomic approaches are providing important insights into nematodes parasitic biology as well as parasite-host interactions. In this scenario, small non-coding RNAs are a promising category of molecules with biological functions 16,33 . Among these, helminths secreted miRNAs could be fundamental in shaping host-parasite interactions and are potential targets for diagnosis and therapy 14 . Some of them are selectively packaged into extracellular vesicles. The interest in EVs is growing due to the increasing studies demonstrating their crucial involvement in intercellular communication, modulation of immune responses, pathology, and their fundamental regulatory role in host-parasite interface during helminthic infections 14 . Parasitic helminths are of particular interest, because they usually determine chronic but rarely lethal infections by manipulating host immune system, thus allowing the www.nature.com/scientificreports/ survival in an adverse niche. Despite their socioeconomic and medical impact, research on this field is still on its infancy, and most of the efforts to characterize nematode miRNAs focused on whole-worm extracts and not on miRNA populations secreted by parasites.
Here, we provide the first miRNAs catalogue from the zoonotic species A. pegreffii with 156 predicted matures, a widespread parasitic nematode of medical and economic concerns. Moreover, we confirm the ability of anisakids to releases EVs during its infective stage (L3). Size of A. pegreffii EVs are in agreement with the only available evidence about the genus Anisakis sp. 22  . Comparisons with material retrievable from public repositories showed 43% (n = 68) of A. pegreffii miRNA catalogue with high level of homology with the phylogenetically related species as A. suum and in around 30% (n = 47) a partial conservation of sequence also with other helminths or organisms was observed. The remaining 25% of miRNAs (n = 41) were observed only in this parasitic organism based on current sequence data, and did not have homologues in any species, being putatively classified as "bona fide" species-specific for A. pegreffii. It is possible that additional non-conserved miRNAs could be present in A. pegreffii and these could be identified by analysing additional life stages and/or with the support of a reliable species-specific genome. In fact, the A. pegreffii reads from RNA-seq were, of necessity, mapped to the A. simplex genome, therefore additional A. pegreffii-specific sequences may have been undetected in the proposed miRNAs list. Similar issues were observed in previous studies we performed for the entire repertoire of transcripts of A. simplex sensu stricto and A. pegreffii, showing ranges of only 80-85% of genome mapping 38,39 , even for specimens belonging to the same species, thus suggesting the low annotation reliability of the available A. simplex genome.
Among the conserved miRNAs, one of the two common miRNAs of nematodes with crucial role in development (lin-4 and let-7) was not observed here: let-7 was not included in the final version of the catalogue according  www.nature.com/scientificreports/ to its very low abundance and its amount of variation compared to known let-7 sequences. A similar evidence was reported in H. contortus 40 . The scientific community is still questioning whether miRNAs from body fluids play biological functions, as the concentration of extracellular miRNAs may be too low to exert in vivo effect 41 .
Half of miRNAs enriched in the extracellular vesicles fraction were novel. Since previous works used the number of sequence reads of a particular miRNA from deep-sequencing as an indication of molecular abundance 40,42 , some miRNAs found at high abundance in infective third-stage larvae were abundant also in the extracellular vesicles fraction (miR-100a-5p, miR-1-3p, miR-71-5p, si-miR-9-5p). In particular, ape-miR-100a-5p is the only miRNAs among the four highly abundant that is upregulated in EVs (FC = 2.49, FDR = 0.00023) and such large presence may indicate a biological role. Another abundant miRNA is ape-lin-4, found both in larvae and EVs. It has been suggested that a dysregulation of miR-100-5p may be related with several types of human cancers [43][44][45]. Similarly, the potential targets of ape-lin-4 are genes involved in cellular proliferation, tumor suppression and induction of apoptosis.
Other miRNAs found differentially expressed and enriched in EVs were mostly novel and not very abundant, but their validity confirmed by Stem and Loop RT-PCR and their selective package into EVs may suggest a potential role in host-parasite interplay. Previous studies identified miR-71 and miR-100c as common markers in the sera from hosts infected with B. malayi, Dirofilaria immitis and Loa loa 34 .
Studies on exosomal miRNAs from human parasitic nematodes have been conducted only on adults and L3 of Brugia 34,46 and Ascaris 37 . Promising evidences were obtained with animal models and a role of exosomal miRNAs in targeting host genes associated with immunity and inflammation has been identified for the intestinal nematode of rodents Heligmosomoides 47 . In fact, administration of nematode exosomes to mice suppresses Type 2 innate responses and eosinophilia induced by the allergenic fungus Alternaria, by silencing Il33r and Dusp1. Table 4. List of ten selected Anisakis pegreffii abundant miRNAs in larvae and in extracellular vesicles, together with putative orthologues miRNAs from other parasitic helminths with conserved seed region, and human miRNAs putative orthologues. The last two columns include Anisakis pegreffii miRNAs predictive targets in human genome according to miRDB and their ID according to NCBI database. Asterisks indicate a miRNAs enriched in exosomes, according to literature. (Underlined miRNAs are abundant both in L3 and in EVs list). Target-putative role  GeneID NCBI   miR-100a-5p asu-miR-100b-5p bma-miR-100a bma-miR-100b bma-miR-100c bma-miR-100d hpo-miR-100-5p* hsa-miR-100-5p
Here, the experimental settings for EVs release by L3 were selected to mimic human body/host conditions (body temperature 37 °C). We would like to point out that humans are accidental hosts, and therefore interactions taking place between Anisakis and humans are not the result of natural evolutionary processes, such as co-evolution and/or co-adaptation. Nevertheless, extracellular vesicles miRNAs released by Anisakis L3 may be involved in pathogenic condition and clinical outcomes, and the analyses of putative miRNAs targets may shed light in their potential pathogenic effect.
Most of the identified miRNAs were common in the two classes of material analysed, namely the infective larval stage and its released extracellular vesicles, as expected. However, a group of 12 miRNAs was enriched in the latter and gene target analyses using human genome as query revealed several predicted gene targets of high interest. These are involved in crucial processes during infections, as cellular proliferation and/or differentiation during the shift from innate to adaptive immune response, apoptosis and inflammation, commonly associated with the functions of extracellular vesicles and exosomes 14 .
In conclusion, the present study provides the first evidence of EVs miRNAs content released by A. pegreffii infective larvae: altogether, the results suggest larval and EVs associated miRNAs may play an important role in the host-parasite interplay, given the wide conservation of sequences among parasitic helminths.

Methods
Fish and parasites collection. Third stage larvae (= L3) of parasitic nematodes belonging to Anisakis genus were collected from fishes purchased at market in 2019-2020. No live fishes have been involved in the study. In brief, visceral body cavity of 10 specimens of Merluccius merluccius and 50 specimens of Engraulis encrasicolus from FAO area 37 (Mediterranean Sea) were inspected and a total of 50 nematodes were selected to perform experiments of small-RNA isolation from L3 and from their released extracellular vesicles (= EVs), collected after incubations procedures described later. The vitality of L3s was evaluated based on their sponta- www.nature.com/scientificreports/ neous movements, after the encapsulation removal. All the L3 selected for the study were mixed from different hosts, in order to have homogeneous samples and avoid any host-related batch effects. Moreover, L3 were identified at species level as belonging to A. pegreffii, using the molecular diagnostic key based on ITS PCR-RFLP procedure 51 . All L3 were washed repeatedly with filtered PBS before molecular procedures.
Size and concentration analysis of EVs. The size distribution and particle concentration of the fraction recovered after the EVs enrichment procedure were measured using Nanoparticle Tracking Analysis. NTA was carried out using a Nanosight NS300 (Malvern Panalytical). Five measurements were performed with 60 s duration of each measurement and the data was analysed using NTA software version 3.4.
RNA extraction, library preparation and RNA-sequencing. Three biological samples of L3 (5 pooled larve) and EVs (5 pooled larve) were used. Total small-RNA fraction from L3 were isolated using miRNeasy tissue kit (Qiagen, Hilden, Germany), according to manufacture instructions. EVs were obtained after incubating pool of L3 in 24 multiwell plates with filtered 1 ml of RPMI with 1% pen-strep for 24 h at 37 °C and 5%CO 2 . Surnatants were collected and Exoquick (System Bio) was used according to manufacturer's instructions to obtain the exosome-enriched fractions (resuspended in 100ul of 20um filtered PBS). From these fractions, small-RNAs were obtained using miRNeasy Serum/Plasma kit (Qiagen, Hilden, Germany Reads aligned to AS14 were also mapped to a collection of Anisakis spp. rRNA and ncRNAs (excluding hairpins and predicted miRNAs) retrieved from the available repositories GenBank, RNA Central and Rfam platform 56 . Reads unaligned to ribosomal RNAs were analysed for their size distribution as reported in Fig. 2 and then mapped (-n 0 -l 18 -a -best -strata -e 80 -norc) to two different lists. The first was composed by 153 putative miRNA precursors (see below) plus other Anisakis ncRNAs. Reads aligning to this list were used for the correlation analysis. The second list included 206 putative mature Anisakis miRNAs (compiled as described below); reads mapping to these miRNAs were used for the differential expression analysis and the heatmap construction.

Prediction of Anisakis pegreffii miRNAs. In the absense of any previous information on miRNAs from
Anisakis species, two independent and complimentary approaches were used to compile a list of putative A. pegreffii miRNAs to be used for mapping. In a first approach, precursors and mature miRNAs from Ascaris suum were retrieved from miRBase (release 22) and used to search potential orthologues in the genome of A. simplex by using the BLAST tool implemented in WormBase Parasite. Ascaris suum was chosen because represented the evolutionary closer parasitic nematode for which miRNA sequence information was available in public repositories, and A. simplex is the only species from the genus Anisakis with a sequenced genome. Hairpins were included according to inclusive parameters, when mature of the two species showed ≥ 85% identity on ≥ 80% length and max 2 mm in the seed (total mm max 6). This way a first list ("dataset1") of Anisakis putative mature and precursor miRNAs was obtained. In a second approach, reads from our samples mapping to the A. simplex genome, and subtracted of those representing rRNAs and other ncRNAs (i.e., tRNAs, snoRNAs, snRNAs), were used to search the A. simplex genome by the miRNA prediction software miRDeep* 57 . This way a second list ("dataset2") was obtained. The two datasets were compared in order to create an inclusive non-redundant inventory of putative hairpins (153) and mature (206) Anisakis miRNAs. Such lists were used for the different following analyses as described. Secondary structure predictions and minimal free energy calculations were performed using RNAfold 58 .
Quantification and differential expression of miRNAs. Reads from A. pegreffii small-RNA libraries were mapped to the final list including 206 putative mature miRNAs and then used for the differential expression analysis. Reads with multiple highest score mappings were discarded. Expression values were calculated as count per millions (CPM, that is the number of reads mapping on a feature divided by the total number of mapped reads and multiplied by one million) and used for sample clustering. Reads mapping to miRNAs and hairpins included in the two datasets were assigned to mature miRNAs. Differential expression analysis of mature miRNAs with 1 CPM in at least two samples of each triplicate was performed using glmFIT and glmLRT functions provided by the edgeR software package 59,60 . Log2 Fold change (FC) and false discovery rates (FDR) were calculated to provide statistical validation (Supplementary Table 2) and inclusive parameters with moderate thresholds were considered for further analyses on DEGs miRNAs (FC > 1 and FDR < 0.05).  61 using as template small-RNA L3 and enriched exosomal fractions of A. pegreffii. First-strand cDNA was generated using the SuperScript II Reverse Transcriptase (Invitrogen) according to manufacturer's instruction and specific stem-loop primers (0.1 μM). Real time PCR amplifications were performed in a final volume of 20 μl including SYBR Green Master Mix (Applied Biosystem), specific forward and universal reverse primers (0.1 μM each) and 2 μl of cDNA. Amplification was as follows: initial holding stage of 2 min at 50 °C and 2 min at 95 °C followed by 40 cycles (30 s. 95 °C, 1 min 60 °C). All RT-qPCR reactions were performed in biological and technical triplicates. The subset of 10 miRNAs was selected according to their total abundance and their enrichment in EVs (ape-novel-miR-19, ape-miR-7-5p, ape-novel-miR-65, ape-novel-miR-27, ape-novel-miR-131, ape-miR-72-5p, ape-novel-miR-184, ape-miR-100a-5p, ape-lin-4-5p, ape-miR-1-3p). The sequence of all primers employed is provided in Supplementary Tables 4a and 4b.
Target analysis and seed conservation. Prediction of putative A. pegreffii miRNAs gene targets was inferred using miRDB implementing MiRTarget bioinformatic tool 62 , an online database for miRNA target prediction and functional annotations, using human genome as query and indicating the related orthologue human miRNA. The top20 most abundant miRNAs and on those enriched in the exosomal fraction were considered for this analysis and only the targets with highest score and potential role in immune or inflammatory response were retained.
With the aim to explore conservation of seed region across parasitic helminthic lineages, we explored the amount of variation in A. pegreffii miRNAs in comparison to those available from other parasitic helminths, according to the supplementary data retrieved 21 , retaining only those with complete seed identity or maximum 1 mismatch. Complete list is available in Supplementary Table 3.

Data availability
The small RNA-Seq datasets generated and analysed in the current study have been deposited in NCBI's Gene Expression Omnibus and are accessible as Bioproject PRJNA786753. Other data generated during this study have been included as Supplementary Information.